function model= kTuning(model, epsilon)

k = model.k;
S = model.Sigma;
w = model.weight;
temp = zeros(k,k);
for m = 1:k
    for n = 1:k
        temp(m,n) = abs(S(:,:,m) - S(:,:,n))/(S(:,:,m)+S(:,:,n));
    end
end

temp = temp + 10000 * eye(k); 

if min(temp(:)) <epsilon
    [pair,~] = find(temp == min(min(temp)));
    i = pair(1);
    j = pair(2);
    w_new = w(i) + w(j);
    S_new = (w(i)*S(:,:,i) + w(j)*S(:,:,j))/(w(i)+w(j));
    % delete original and add the newly generated one to the end
    w([i,j]) = [];
    weight = [w,w_new];
    S(:,:,[i,j]) = [];
    Sigma = cat(3,S,S_new);
    
    model.k = k-1;
    model.Sigma = Sigma;
    model.weight = weight;
end
